mvgam case study 2: multivariate models
Nicholas Clark (n.clark@uq.edu.au)
In this example we will examine multivariate forecasting models using
mvgam, which fits dynamic GAMs (DGAMs) using MCMC sampling
(note that either JAGS or Stan is required;
installation links are found here and here). First a
simulation experiment to determine whether mvgam‘s
inclusion of complexity penalisation works by reducing the number of
un-needed dynamic factors. In any factor model, choosing the appropriate
number of factors M can be difficult. The approach used by
mvgam when sampling with JAGS is to estimate a
penalty for each factor that squeezes the factor’s variance toward zero,
effectively forcing the factor to evolve as a flat white noise process.
By allowing each factor’s penalty to be estimated in an exponentially
increasing manner (following Welty, Leah J., et al. Bayesian distributed
lag models: estimating effects of particulate matter air pollution on
daily mortality Biometrics 65.1 (2009): 282-291), we hope that we can
guard against specifying too large a M. Note that when
sampling with Stan, we capitalise on the superior sampling
and exploration of Hamiltonian Monte Carlo to choose the number of
factors by placing independent normal priors on factor standard
deviations. Begin by simulating 6 series that evolve with a
shared seasonal pattern and that depend on 2 latent random
walk factors. Each series is 100 time steps long, with a
seasonal frequency of 12. We give the trend moderate
importance by setting trend_rel = 0.6 and we allow each
series’ observation process to be drawn from slightly different Poisson
distributions
set.seed(1111)
library(mvgam)## Loading required package: mgcv## Loading required package: nlme## This is mgcv 1.8-33. For overview type 'help("mgcv-package")'.## Loading required package: paralleldat <- sim_mvgam(T = 100, n_series = 6, n_lv = 2,
family = "poisson", mu_obs = runif(8,
4, 6), trend_rel = 0.6, train_prop = 0.85)## Warning: replacing previous import 'lifecycle::last_warnings' by
## 'rlang::last_warnings' when loading 'tibble'## Warning: replacing previous import 'lifecycle::last_warnings' by
## 'rlang::last_warnings' when loading 'pillar'## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zooHave a look at the series
par(mfrow = c(3, 2))
for (i in 1:6) {
plot(dat$data_train$y[which(as.numeric(dat$data_train$series) ==
i)], type = "l", ylab = paste("Series",
i), xlab = "Time", bty = "L")
box(bty = "L", lwd = 2)
}
par(mfrow = c(1, 1))Clearly there are some correlations in the trends for these series.
But how does a dynamic factor process allow us to potentially capture
these dependencies? The below example demonstrates how. Essentially, a
dynamic factor is an unobserved (latent) random process that
induces correlations between time series via a set of factor loadings
(\(\beta\)) while exercising dimension
reduction. The loadings represent constant associations between the
observed time series and the dynamic factor, but each series can still
deviate from the factor through its error process and its associations
with other factors (if we estimate >1 latent factor in
our model).
A challenge with any factor model is the need to determine the number
of factors M. Setting M too small prevents
temporal dependencies from being adequately modelled, leading to poor
convergence and difficulty estimating smooth parameters. By contrast,
setting M too large leads to unnecessary computation.
mvgam approaches this problem by formulating a prior
distribution that enforces exponentially increasing penalties on the
factor variances to allow any un-needed factors to evolve as flat lines.
Now let’s fit a well-specified model for our simulated series in which
we estimate random intercepts, a shared seasonal cyclic smooth and
2 latent dynamic factors
mod1 <- mvgam(data_train = dat$data_train,
data_test = dat$data_test, formula = y ~
s(series, bs = "re") + s(season,
bs = c("cc"), k = 8), knots = list(season = c(0.5,
12.5)), use_lv = TRUE, n_lv = 2,
family = "poisson", trend_model = "RW",
use_stan = TRUE, chains = 4, burnin = 500)## Info: Found int division at 'string', line 22, column 31 to column 48:
## n_lv * (n_lv - 1) / 2
## Values will be rounded towards zero. If rounding is not desired you can write
## the division as
## n_lv * (n_lv - 1) / 2.0Look at a few plots. The estimated smooth function
plot_mvgam_smooth(object = mod1, series = 1,
smooth = "season")
And the true seasonal function in the simulation
plot(dat$global_seasonality[1:12], type = "l",
bty = "L", ylab = "True function", xlab = "season")
Check whether each factor was retained using the
plot_mvgam_factors function. Here, each factor is tested
against a null hypothesis of white noise by calculating the sum of the
factor’s 1st derivatives. A factor that has a larger contribution to the
series’ latent trends will have a larger sum, both because that factor’s
absolute magnitudes will be larger (due to the weaker penalty on the
factor’s precision) and because the factor will move around more. By
normalising these estimated first derivative sums, it should be apparent
whether some factors have been dropped from the model. Here we see that
each factor is contributing to the series’ latent trends, and the plots
show that neither has been forced to evolve as white noise
plot_mvgam_factors(mod1)
Now we fit the same model but assume that we no nothing about how
many factors to use, so we specify the maximum allowed (the total number
of series; 6). Note that this model is computationally more
expensive so it will take longer to fit
mod2 <- mvgam(data_train = dat$data_train,
data_test = dat$data_test, formula = y ~
s(series, bs = "re") + s(season,
bs = c("cc"), k = 8), knots = list(season = c(0.5,
12.5)), use_lv = TRUE, n_lv = 6,
family = "poisson", use_stan = TRUE,
trend_model = "RW", chains = 4, burnin = 500)## If rounding is intended please use the integer division operator %/%.
## Info: Found int division at 'string', line 22, column 31 to column 48:
## n_lv * (n_lv - 1) / 2
## Values will be rounded towards zero. If rounding is not desired you can write
## the division as
## n_lv * (n_lv - 1) / 2.0Use the same plots as model 1 to see if this model has also fit the data well
plot_mvgam_smooth(object = mod2, series = 1,
smooth = "season")
plot(dat$global_seasonality[1:12], type = "l",
bty = "L", ylab = "True function", xlab = "season")
Examining the factor contributions gives us some insight into whether
we set n_lv larger than we perhaps needed to (with some of
the factors clearly evolving as unconstrained, zero-centred random
walks). These contributions can be interpreted similarly to ordination
axes when deciding how many latent variables to specify
plot_mvgam_factors(mod2)
The very weak contributions by some of the factors are a result of
the penalisation, which will become more important as the dimensionality
of the data grows. Now onto an empirical example. Here we will access
monthly search volume data from Google Trends, focusing on
relative importances of search terms related to tick paralysis in
Queensland, Australia
library(tidyr)
if (!require(gtrendsR)) {
install.packages("gtrendsR")
}
terms = c("tick bite", "tick paralysis",
"dog tick", "paralysis tick dog")
trends <- gtrendsR::gtrends(terms, geo = "AU-QLD",
time = "all", onlyInterest = T)Google Trends modified their algorithm for extracting
search volume data in 2012, so we filter the series to only include
observations after this point in time
gtest <- trends$interest_over_time %>% tidyr::spread(keyword,
hits) %>% dplyr::select(-geo, -time,
-gprop, -category) %>% dplyr::mutate(date = lubridate::ymd(date)) %>%
dplyr::mutate(year = lubridate::year(date)) %>%
dplyr::filter(year > 2012) %>% dplyr::select(-year)Convert to an xts object and then to the required
mvgam format, holding out the final 10% of observations as
the test data
series <- xts::xts(x = gtest[, -1], order.by = gtest$date)
trends_data <- series_to_mvgam(series, freq = 12,
train_prop = 0.9)Plot the series to see how similar their seasonal shapes are over time
plot(series, legend.loc = "topleft")
Now we will fit an mvgam model with shared seasonality
and random intercepts per series. Our first attempt will ignore any
temporal component in the residuals so that we can identidy which GAM
predictor combination gives us the best fit, prior to investigating how
to deal with any remaining autocorrelation. We assume a Poisson
observation model for the response. Also note that any smooths using the
random effects basis (s(series, bs = "re") below) are
automatically re-parameterised to use the non-centred
parameterisation that is necessary to help avoid common posterior
degeneracies in hierarchical models. This parameterisation tends to
work better for most ecological problems where the data for each group /
context are not highly informative, but it is still probably worth
investigating whether a centred or even a mix of centred / non-centred
will give better computational performance. We suppress the global
intercept as it is not needed and will lead to identifiability issues
when estimating the series-specific random intercepts
trends_mod1 <- mvgam(data_train = trends_data$data_train,
data_test = trends_data$data_test, formula = y ~
s(series, bs = "re") + s(season,
k = 8, m = 2, bs = "cc") - 1,
knots = list(season = c(0.5, 12.5)),
trend_model = "None", family = "poisson",
use_stan = TRUE, chains = 4, burnin = 1000)Given that these series could potentially be following a hierarchical seasonality, we will also trial a slghtly more complex model with an extra smoothing term per series that allows its seasonal curve to deviate from the global seasonal smooth. Ignore the warning about repeated smooths, as this is not an issue for estimation.
trends_mod2 <- mvgam(data_train = trends_data$data_train,
data_test = trends_data$data_test, formula = y ~
s(season, k = 8, m = 2, bs = "cc") +
s(season, series, k = 5, bs = "fs",
m = 1), knots = list(season = c(0.5,
12.5)), trend_model = "None", family = "poisson",
use_stan = TRUE, chains = 4, burnin = 1000)How can we compare these models to ensure we choose one that performs
well and provides useful inferences? Beyond posterior retrodictive and
predictive checks, we can take advantage of the fact that
mvgam fits an mgcv model to provide all the
necessary penalty matrices, as well as to identify good initial values
for smoothing parameters. Because we did not modify this model by adding
a trend component (the only modification is that we estimated
series-specific overdispersion parameters), we can still employ the
usual mgcv model comparison routines
anova(trends_mod1$mgcv_model, trends_mod2$mgcv_model,
test = "LRT")AIC(trends_mod1$mgcv_model, trends_mod2$mgcv_model)summary(trends_mod1$mgcv_model)##
## Family: poisson
## Link function: log
##
## Formula:
## y ~ s(series, bs = "re") + s(season, k = 8, m = 2, bs = "cc") -
## 1
##
## Approximate significance of smooth terms:
## edf Ref.df Chi.sq p-value
## s(series) 4.000 4 67947 <2e-16 ***
## s(season) 5.646 6 4464129 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.795 Deviance explained = 81.4%
## -REML = 1392.2 Scale est. = 1 n = 412summary(trends_mod2$mgcv_model)##
## Family: poisson
## Link function: log
##
## Formula:
## y ~ s(season, k = 8, m = 2, bs = "cc") + s(season, series, k = 5,
## bs = "fs", m = 1)
##
## Parametric coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.6755 0.4403 6.077 1.23e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df Chi.sq p-value
## s(season) 5.507 6 177.7 <2e-16 ***
## s(season,series) 14.165 19 1851.9 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.831 Deviance explained = 85%
## -REML = 1322 Scale est. = 1 n = 412Model 2 seems to fit better so far, suggesting that hierarchical
seasonality gives better performance for these series. But a problem
with both of the above models is that their forecast uncertainties will
not increase into the future, which is not how time series forecasts
should behave. Here we fit Model 2 again but specifying a time series
model for the latent trends. We assume the dynamic trends can be
represented using latent factors that each follow a RW process, and we
will rely on the exponential penalties to help regularise any un-needed
factors by setting n_lv = 4
trends_mod3 <- mvgam(data_train = trends_data$data_train,
data_test = trends_data$data_test, formula = y ~
s(season, k = 8, m = 2, bs = "cc") +
s(season, series, k = 5, bs = "fs",
m = 1), knots = list(season = c(0.5,
12.5)), trend_model = "RW", use_lv = TRUE,
n_lv = 4, family = "poisson", use_stan = TRUE,
chains = 4, burnin = 1000)## If rounding is intended please use the integer division operator %/%.
## Info: Found int division at 'string', line 24, column 31 to column 48:
## n_lv * (n_lv - 1) / 2
## Values will be rounded towards zero. If rounding is not desired you can write
## the division as
## n_lv * (n_lv - 1) / 2.0Have a look at the returned Stan model file to see how
the dynamic factors are incorporated
trends_mod3$model_file## [1] ""
## [2] "//Stan model code generated by package mvgam"
## [3] "data {"
## [4] "int idx1[16]; // discontiguous index values"
## [5] "int idx2[4]; // discontiguous index values"
## [6] "int<lower=0> total_obs; // total number of observations"
## [7] "int<lower=0> n; // number of timepoints per series"
## [8] "int<lower=0> n_lv; // number of dynamic factors"
## [9] "int<lower=0> n_sp; // number of smoothing parameters"
## [10] "int<lower=0> n_series; // number of series"
## [11] "int<lower=0> num_basis; // total number of basis coefficients"
## [12] "real p_taus[1]; // prior precisions for parametric coefficients"
## [13] "real p_coefs[1]; // prior locations for parametric coefficients"
## [14] "vector[num_basis] zero; // prior locations for basis coefficients"
## [15] "matrix[num_basis, total_obs] X; // transposed mgcv GAM design matrix"
## [16] "int<lower=0> ytimes[n, n_series]; // time-ordered matrix (which col in X belongs to each [time, series] observation?)"
## [17] "matrix[6,6] S1; // mgcv smooth penalty matrix S1"
## [18] "int<lower=0, upper=1> y_observed[n, n_series]; // indices of missing vs observed"
## [19] "int<lower=-1> y[n, n_series]; // time-ordered observations, with -1 indicating missing"
## [20] "}"
## [21] ""
## [22] "transformed data{"
## [23] "// Number of non-zero lower triangular factor loadings"
## [24] "// Ensures identifiability of the model - no rotation of factors"
## [25] "int<lower=1> M;"
## [26] "M = n_lv * (n_series - n_lv) + n_lv * (n_lv - 1) / 2 + n_lv;"
## [27] "}"
## [28] ""
## [29] "parameters {"
## [30] "// raw basis coefficients"
## [31] "row_vector<lower=-30,upper=30>[num_basis] b_raw;"
## [32] ""
## [33] "// dynamic factors"
## [34] "matrix[n, n_lv] LV;"
## [35] ""
## [36] "// dynamic factor lower triangle loading coefficients"
## [37] "vector[M] L;"
## [38] ""
## [39] "// smoothing parameters"
## [40] "vector<lower=0.0005>[n_sp] lambda;"
## [41] "}"
## [42] ""
## [43] "transformed parameters {"
## [44] "// GAM contribution to expectations (log scale)"
## [45] "vector[total_obs] eta;"
## [46] ""
## [47] "// trends and dynamic factor loading matrix"
## [48] "matrix[n, n_series] trend;"
## [49] "matrix[n_series, n_lv] lv_coefs;"
## [50] ""
## [51] "// basis coefficients"
## [52] "row_vector[num_basis] b;"
## [53] ""
## [54] "for (i in 1:num_basis) {"
## [55] "b[i] = b_raw[i];"
## [56] "}"
## [57] "// constraints allow identifiability of loadings"
## [58] "for (i in 1:(n_lv - 1)) {"
## [59] "for (j in (i + 1):(n_lv)){"
## [60] "lv_coefs[i, j] = 0;"
## [61] "}"
## [62] "}"
## [63] "{"
## [64] "int index;"
## [65] "index = 0;"
## [66] "for (j in 1:n_lv) {"
## [67] "for (i in j:n_series) {"
## [68] "index = index + 1;"
## [69] "lv_coefs[i, j] = L[index];"
## [70] "}"
## [71] "}"
## [72] "}"
## [73] ""
## [74] "// derived latent trends"
## [75] "for (i in 1:n){;"
## [76] "for (s in 1:n_series){"
## [77] "trend[i, s] = dot_product(lv_coefs[s,], LV[i,]);"
## [78] "}"
## [79] "}"
## [80] ""
## [81] "eta = to_vector(b * X);"
## [82] "}"
## [83] ""
## [84] "model {"
## [85] "// parametric effect priors (regularised for identifiability)"
## [86] "for (i in 1:1) {"
## [87] "b_raw[i] ~ normal(p_coefs[i], 1 / p_taus[i]);"
## [88] "}"
## [89] ""
## [90] "// prior for s(season)..."
## [91] "b_raw[2:7] ~ multi_normal_prec(zero[2:7],S1[1:6,1:6] * lambda[1]);"
## [92] ""
## [93] "// prior for s(season,series)..."
## [94] "for (i in idx1) { b[i] ~ normal(0, lambda[2]); }"
## [95] ""
## [96] "for (i in idx2) { b[i] ~ normal(0, lambda[3]); }"
## [97] ""
## [98] "// priors for smoothing parameters"
## [99] "lambda ~ exponential(0.05);"
## [100] ""
## [101] "// priors for dynamic factor loading coefficients"
## [102] "L ~ double_exponential(0, 1);"
## [103] ""
## [104] "// dynamic factor estimates"
## [105] "for (j in 1:n_lv) {"
## [106] "LV[1, j] ~ normal(0, 1);"
## [107] "}"
## [108] ""
## [109] "for (j in 1:n_lv) {"
## [110] "LV[2:n, j] ~ normal(LV[1:(n - 1), j], 1);"
## [111] "}"
## [112] ""
## [113] "// likelihood functions"
## [114] "for (i in 1:n) {"
## [115] "for (s in 1:n_series) {"
## [116] "if (y_observed[i, s])"
## [117] "y[i, s] ~ poisson_log(eta[ytimes[i, s]] + trend[i, s]);"
## [118] "}"
## [119] "}"
## [120] "}"
## [121] ""
## [122] "generated quantities {"
## [123] "vector[n_sp] rho;"
## [124] "vector[n_lv] penalty;"
## [125] "matrix[n, n_series] ypred;"
## [126] "rho = log(lambda);"
## [127] "penalty = rep_vector(1.0, n_lv);"
## [128] ""
## [129] "// posterior predictions"
## [130] "for(i in 1:n){"
## [131] "for(s in 1:n_series){"
## [132] "ypred[i, s] = poisson_log_rng(eta[ytimes[i, s]] + trend[i, s]);"
## [133] "}"
## [134] "}"
## [135] "}"
## [136] ""Inspection of the dynamic factors and their relative contributions indicates that the first factor is by far the most important
plot_mvgam_factors(trends_mod3)
Model 3 (with the dynamic trend) should provide far superior forecasts than relying only on the estimated smooths. Inspect the model summary
summary(trends_mod3)## GAM formula:## y ~ s(season, k = 8, m = 2, bs = "cc") + s(season, series, k = 5,
## bs = "fs", m = 1)## ## Family:## Poisson## ## Link function:## log## ## Trend model:## RW## ## N latent factors:## 4## ## N series:## 4## ## N observations:## 412## ## Status:## Fitted using rstan::stan()## ## GAM smooth term estimated degrees of freedom:## edf df
## s(season) 15.34 6
## s(season,series) 241.31 5## ## GAM coefficient (beta) estimates:## 2.5% 50% 97.5% Rhat n.eff
## (Intercept) 2.015451599 2.704062661 3.357046358 1 3827
## s(season).1 -0.725954600 -0.558499214 -0.391935085 1 4428
## s(season).2 -0.745968815 -0.530358930 -0.331730807 1 4622
## s(season).3 -0.454880908 -0.243435772 -0.034444696 1 4279
## s(season).4 -0.126528126 0.096704318 0.298432717 1 3892
## s(season).5 0.594950983 0.795557653 0.996021247 1 4303
## s(season).6 0.415311823 0.570308255 0.727358055 1 4434
## s(season,series).1 -0.025902484 0.112598184 0.270454054 1 4321
## s(season,series).2 -0.037057937 0.116948837 0.282879379 1 4343
## s(season,series).3 -0.235609594 -0.088417101 0.060802425 1 4144
## s(season,series).4 0.075293472 0.146413231 0.226387481 1 4321
## s(season,series).5 -0.246916278 0.481381022 1.224551401 1 2620
## s(season,series).6 -0.349067798 -0.141097228 0.039786395 1 4101
## s(season,series).7 -0.146865922 0.049291945 0.253029573 1 4909
## s(season,series).8 -0.197040972 -0.032679354 0.134856058 1 4284
## s(season,series).9 -0.090271168 0.001148285 0.094937042 1 4422
## s(season,series).10 -1.770867701 -0.993878954 -0.223043479 1 2947
## s(season,series).11 -0.086619972 0.065147423 0.223972260 1 4423
## s(season,series).12 -0.273598256 -0.104589717 0.054089139 1 4330
## s(season,series).13 -0.004182408 0.135707321 0.295333196 1 3979
## s(season,series).14 0.034888627 0.109534198 0.191690222 1 4289
## s(season,series).15 -0.929026883 -0.198378160 0.527048316 1 3431
## s(season,series).16 -0.308198606 -0.147532763 -0.002552626 1 4508
## s(season,series).17 -0.170926173 -0.010173041 0.139672936 1 4541
## s(season,series).18 -0.153108141 -0.006980581 0.139024325 1 4691
## s(season,series).19 -0.063331002 0.006514286 0.080009901 1 4444
## s(season,series).20 0.122248054 0.840516926 1.539518362 1 3712## ## GAM smoothing parameter (rho) estimates:## 2.5% 50% 97.5% Rhat n.eff
## s(season) 2.3903058 3.610948727 4.490441 1 4569
## s(season,series) -2.5442110 -2.075032898 -1.560145 1 3100
## s(season,series)2 -0.7368928 0.006311793 1.167338 1 3516## ## [1] "n_eff / iter looks reasonable for all parameters"
## [1] "Rhat looks reasonable for all parameters"
## [1] "0 of 4000 iterations ended with a divergence (0%)"
## [1] "0 of 4000 iterations saturated the maximum tree depth of 12 (0%)"
## [1] "E-FMI indicated no pathological behavior"Look at Dunn-Smyth residuals for some series from this preferred model to ensure that our dynamic factor process has captured most of the temporal dependencies in the observations
plot_mvgam_resids(trends_mod3, series = 1)
plot_mvgam_resids(trends_mod3, series = 2)
plot_mvgam_resids(trends_mod3, series = 3)
plot_mvgam_resids(trends_mod3, series = 4)
Perform posterior predictive checks to see if the model is able to
simulate data that looks realistic and unbiased by examining simulated
kernel densities for posterior predictions (yhat) compared
to the density of the observations (y). This will be
particularly useful for examining whether the Negative Binomial
observation model is able to produce realistic looking simulations for
each individual series.
ppc(trends_mod3, series = 1, type = "hist")
ppc(trends_mod3, series = 2, type = "hist")
ppc(trends_mod3, series = 3, type = "hist")
ppc(trends_mod3, series = 4, type = "hist")
Look at traceplots for the smoothing parameters
(rho)
plot_mvgam_trace(object = trends_mod3, param = "rho")
Plot posterior predictive distributions for the training and testing periods for each series
plot_mvgam_fc(object = trends_mod3, series = 1,
data_test = trends_data$data_test)## Out of sample DRPS:## [1] 51.25285## 
plot_mvgam_fc(object = trends_mod3, series = 2,
data_test = trends_data$data_test)## Out of sample DRPS:## [1] 13.42173## 
plot_mvgam_fc(object = trends_mod3, series = 3,
data_test = trends_data$data_test)## Out of sample DRPS:## [1] 43.18961## 
plot_mvgam_fc(object = trends_mod3, series = 4,
data_test = trends_data$data_test)## Out of sample DRPS:## [1] 33.33799## 
Plot posterior distributions for the latent trend estimates, again for the training and testing periods
plot_mvgam_trend(object = trends_mod3, series = 1,
data_test = trends_data$data_test)
plot_mvgam_trend(object = trends_mod3, series = 2,
data_test = trends_data$data_test)
plot_mvgam_trend(object = trends_mod3, series = 3,
data_test = trends_data$data_test)
plot_mvgam_trend(object = trends_mod3, series = 4,
data_test = trends_data$data_test)
Given that we fit a model with hierarchical seasonality, the seasonal
smooths are able to deviate from one another (though they share the same
wiggliness and all deviate from a common ‘global’ seasonal function).
Here we use the newdata argument to generate predictions
for each of the hierarchical smooth functions (note that the intercept
is still included in these plots so they do not center on zero)
newdat <- data.frame(season = seq(1, 12,
length.out = 100), series = levels(trends_data$data_train$series)[1])
plot_mvgam_smooth(object = trends_mod3, series = 1,
smooth = "season", newdata = newdat)
newdat <- data.frame(season = seq(1, 12,
length.out = 100), series = levels(trends_data$data_train$series)[2])
plot_mvgam_smooth(object = trends_mod3, series = 2,
smooth = "season", newdata = newdat)
newdat <- data.frame(season = seq(1, 12,
length.out = 100), series = levels(trends_data$data_train$series)[3])
plot_mvgam_smooth(object = trends_mod3, series = 3,
smooth = "season", newdata = newdat)
newdat <- data.frame(season = seq(1, 12,
length.out = 100), series = levels(trends_data$data_train$series)[4])
plot_mvgam_smooth(object = trends_mod3, series = 4,
smooth = "season", newdata = newdat)
Plot posterior mean estimates of latent trend correlations. These correlations are more useful than looking at latent factor loadings, for example to inspect ordinations. This is because the orders of the loadings (although constrained for identifiability purposes) can vary from chain to chain
correlations <- lv_correlations(object = trends_mod3)library(ggplot2)## Warning: package 'ggplot2' was built under R version 4.0.5mean_correlations <- correlations$mean_correlations
mean_correlations[upper.tri(mean_correlations)] <- NA
mean_correlations <- data.frame(mean_correlations)
ggplot(mean_correlations %>% tibble::rownames_to_column("series1") %>%
tidyr::pivot_longer(-c(series1), names_to = "series2",
values_to = "Correlation"), aes(x = series1,
y = series2)) + geom_tile(aes(fill = Correlation)) +
scale_fill_gradient2(low = "darkred",
mid = "white", high = "darkblue",
midpoint = 0, breaks = seq(-1, 1,
length.out = 5), limits = c(-1,
1), name = "Trend\ncorrelation") +
labs(x = "", y = "") + theme_dark() +
theme(axis.text.x = element_text(angle = 45,
hjust = 1))
There is certainly some evidence of positive trend correlations for a few of these search terms, which is not surprising given how similar some of them are and how closely linked they should be to interest about tick paralysis in Queensland. Plot some STL decompositions of these series to see if these trends are noticeable in the data
plot(stl(ts(as.vector(series$`tick paralysis`),
frequency = 12), "periodic"))
plot(stl(ts(as.vector(series$`paralysis tick dog`),
frequency = 12), "periodic"))
plot(stl(ts(as.vector(series$`dog tick`),
frequency = 12), "periodic"))
plot(stl(ts(as.vector(series$`tick bite`),
frequency = 12), "periodic"))
Forecast period posterior predictive checks suggest that the model still has room for improvement:
ppc(trends_mod3, series = 1, type = "hist",
data_test = trends_data$data_test)
ppc(trends_mod3, series = 1, type = "mean",
data_test = trends_data$data_test)
ppc(trends_mod3, series = 2, type = "hist",
data_test = trends_data$data_test)
ppc(trends_mod3, series = 2, type = "mean",
data_test = trends_data$data_test)
ppc(trends_mod3, series = 3, type = "hist",
data_test = trends_data$data_test)
ppc(trends_mod3, series = 3, type = "mean",
data_test = trends_data$data_test)
ppc(trends_mod3, series = 4, type = "hist",
data_test = trends_data$data_test)
ppc(trends_mod3, series = 4, type = "mean",
data_test = trends_data$data_test)
Other next steps could involve devising a more goal-specific set of posterior predictive checks (see this paper by Gelman et al and relevant works by Betancourt for examples) and compare out of sample Discrete Rank Probability Scores for this model and other versions for the latent trends (i.e. AR2, AR3, Random Walk)